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ABSTRACT 

We investigate the predictions for the faint-end quasar luminosity function (QLF) and 
its evolution using fully cosmological hydrodynamic simulations which self-consistently 
follow star formation, black hole growth and associated feedback processes. We find 
remarkably good agreement between predicted and observed faint end of the optical 
and X-ray QLFs (the bright end is not accessible in our simulated volumes) at z < 
2. At higher redshifts our simulations tend to overestimate the QLF at the faintest 
luminosities. We show that although the low (high) luminosity ranges of the faint-end 
QLF are dominated by low (high) mass black holes, a wide range of black hole masses 
still contributes to any given luminosity range. This is consistent with the complex 
lightcurves of black holes resulting from the detailed hydrodynamics followed in the 
simulations. Consistent with the results on the QLFs, we find good agreement for the 
evolution of the comoving number density (in optical, soft and hard X-ray bands) of 
AGN for luminosities ^ 10^"^ erg s~^. However, the luminosity density evolution from 
the simulation appears to imply a peak at higher redshift than constrained from hard 
X-ray data (but not in optical). Our predicted excess at the faintest fluxes at z ^ 2 
does not lead to an overestimate to the total X-ray background and its contribution is 
at most a factor of two larger than the unresolved fraction of the 2-8 keV background. 
Even though this could be explained by some yet undetected, perhaps heavily obscured 
faint quasar population, we show that our predictions for the faint sources at high 
redshifts (which are dominated by the low mass black holes) in the simulations are 
likely affected by resolution effects. 

Key words: quasars: general, methods: numerical, black hole physics, galaxies: ac- 
tive, galaxies: nuclei, galaxies: evolution 



1 INTRODUCTION 

In recent years quasars have been used as instrumental tools 
for probing properties of their host galaxies as well as large 
scale structure through cosmic tim e. The existence of black 
holes at the centre of most galaxies l|Kormendy fc Richstonel 
1 19951 ) combined with the correlation between super mas- 
sive b l ack holes and their p a rent g alaxies (Mag orria n ct al.l 
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19981: iFerrarese fc MerrittI l20od: [Ccbhardt et al 



tion 



Tremaine et al.l 20021 : Graham fc Drivern2007il significantly 



strengthen the link between the black hole and the for- 
mation and evolution of galaxies. Although the origins of 
these correlations are not completely understood, recent 
observational and computational studies point to the fun- 
damental role of some form of quasa r feedback for estab- 
lishing them (e.g.lBurkert fc Silkll2q0ll:lGranato et al.ll2004 
ISazonov et al.l [2004 : Springel et al.ll2005a : Churazov et al.l 



One fundamental aspect of the study of quasars is 
the form and evolution of the Quasar Luminosity Func- 
(QLF). Recent surveys, including SDSS (York et al 
and the 2dF QSO Redshift Survey (jLewis et al 
are now providing large samples over sufficient 
ranges that the QLF shape and evolution can 
Also, numerous studies of the 



2000) 
2003), 
redshift 

be investigated in detail. 



QLF have been made , covering the X-ray (jPage et al 
1997]; iMivaii et aD I2001I: iLa Franca et al.ll200l iFiore et al 
2003: lUeda et aLMgOOa: ICowie et all l2003l: 1 




'Barger et al 



20051: iLa Franca et all l2005l: iSilvcrma n et al 
lEbrero et all l2009l : lYencho et all 120091 ) . optical 
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dWolf et aLll2003: ICroom et al.ll2004l: iRichards et al.ll200d ). 



radio llCirasuolo et al l l2005l: (wall et al l l2005h , and IR 
ijMatute et al.l l2006l : iBrown et all l2006l ) bands. Overall 
these studies suggest that the spatial density of quasars 
undergoes a luminosity- dependent evolution, with the 
density of more luminous quasars peaking at higher redshift 
than the less luminous populations. 

Theoretical investigation of t he QLF has been done 
using semi-analyt ical m o dels ( e.g. [Kauffman n fc Haehneltl 
2OOOI: [V olontcr i et al l l2003l: IWv itlic & Loebl l2003l: 
Grana.to et all |2004|: iMalbon et al.l 120071 : iMaruUi et al.1 
20081 : iBonoli et all 120091 ') . Since these models do not self- 



consistently follow black hole growth, the AGN lightcurves 
and luminosity have to be calculated via a specified 
prescription. The predominant method for modeling the 
quasars in this context is to treat quasars as radiating at 
a fixed fraction of their Eddington luminosity for a char- 
acteristic time-scale after a galaxy merger before shutting 
off completely due to feedback efi'ects. The determination 
of the chara cteristic time- s cale varies between methods. 
For example, iHaiman et al. assume quasar radiation 

at the Eddington luminosity for a fixed time-scale of 
2 X 10^ years f or the radio- loud lifetime of the quasars. 
IVolonteri et"al] (|2003D assume the quasar will maintain 
Eddington accretion until it has accreted a total mass 
proportional to the fifth power of the c ircular speed of 
the merged system. IWvithe fc Loebl (l2003t ) adopt a model 
where the quasars radiate at a fixed fraction of the Edding- 
ton luminosity for the dynamical time of the galactic disc, 
at which point the gas has been given more energy than its 
binding energy. These methods have produced promising 
results, but are all based on variants of the simple on-off 
model. 

iHopkins et al] (|2005al lbll3. l2006al lbh took a different ap- 
proach to modeling the QLF by analyzing the light curves 
of quasars in hydrodynamical galaxy merger simulations 
whi ch included black hole growth, accretion an d feedback 
fsee lDi Matteo et al. I I2OO5I : ISpringel et aLll2005bl ). and used 
the results to express the quasar lifetime as a differential 
time a quasar spen ds radiating in a logarithmic luminosity 
bin l|Hopkins et al.. 2005a, b.q"). The quasar lifetimes were fit 
to a Schechter function dependent on both peak and cur- 
rent luminosity. In this way, quasars were modeled using de- 
tailed predictions from hydrodynamic simulations for their 
lightcurves and were shown to radiate at a range of lumi- 
nosities both at and below their peak, rather than being 
restricted to radiating at a primarily constant peak lumi- 
nosity. 

Hopkins et al.'s approach found that using the predicted 
form for quasar lifetime, the faint-end of the QLF could be 
explained by quasars radiating well below their peak lumi- 
nosities, rather than by quasars with low peak luminosities. 
In this case, to match the observational form of the QLF, the 
quasar creation rate must peak at the critical break luminos- 
ity of the QLF, with a very rapid dro p-off for luminosities 
below the break ("Hopk ins et al.ll2005bl ). This work provided 
a fundamentally different explanation for the physical source 
of the faint-end slope and the break luminosity while still 
reproducing the form and evolution of the observed QLF 
l|Hopkins et al.||2006bl ). However, the conclusions were based 
upon data extracted from individual galaxy merger simu- 



lations and have yet to be investigated with cosmological 
hydrodynamical simulations. 

In this paper we analyse fully cosmological hydrody- 
namic simulations which directly include modeling of black 
hole growth, accretion and associated feedback processes (as 
well as the dynamics of dark matter, dissipation, star for- 
mation and stellar feedback) and make predictions for the 
quasar luminosity function and its evolution. The simula- 
tions are currently among the largest, highest- resolution hy- 
drodynamical simulations which include gas hydrodynam- 
ics, and have been shown already to reproduce many as- 
pects of the black hole evolution, such as the mass func- 
tion and accretion rate distribution, and in particular the 
assem bly and evolution of the black hole galaxy correla- 
tions (|Di Matteo et al1l2008l ). In this paper we compare the 
black hole luminosity functions and their evolution from 
the simulations with appropriate observations in various en- 
ergy bands. This is both an important test for assessing the 
value of the simulations and for providing a physical con- 
text within which to interpret the observations and quasar 
evolution in general. 

In Section 2 we describe the numerical modeling for the 
black holes accretion and luminosity (Section 2.1) and the 
simulation parameters used (Section 2.2). In Section 3 we 
present the results for the black hole luminosity function, 
comoving number density evolution, and luminosity density 
evolution, and compare with observational data. In Section 4 
we discuss the implications of our simulation on the hard X- 
ray background, and in Section 5 we summarize and discuss 
our results. 



2 METHOD 

2.1 Numerical simulation 

In this study, we analy se the set of simulations published in 
iDi Matteo et al.l (|2008l ). Here we present a brief summary 
of the si mulation code and the method used. We refer the 
reader to lDi Matteo et al.1 (|2008l ') for all details. 

The code we use is the massively parallel cosmological 
TreePM-SPH code Gadget2 (Springel 2005), with the addi- 
tion of a multi-phase modelling of the ISM, which allows 
treatment of star formation (Springel & Hernquist 2003), 
and black hole accretion and associated feedback processes 
(Springel et al. 2005, Di Matteo et al. 2005). 

Black holes are simulated with coUisionless parti- 
cles that are created in newly emerging and resolved 
groups/galaxies. A friends-of-friends group finder is called 
at regular intervals on the fly (the time intervals are equally 
spaced in log a, with A log a = log 1.25), and employed to 
find groups of particles. Each group that does not already 
contain a black hole is provided with one by turning its 
densest particle into a sink particle with a seed black hole of 
fixed mass, M = 10^/i~^ M0. The black hole particle then 
grows in mass via accretion of surrounding gas according to 

Mbh = 



(|Hovle fc Lvttletonll 19391 : iBondi fc Hovld 



ll944l : lBondiHl952l '). and by merging with other black holes. 

For the simulations used here it is assumed that accre- 
tion is limited to a maximum of 3 times the Eddington rate. 
Note, very few sources accrete at this critical value, as seen 
in Fig. [U 
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Table 1. Numerical Parameters 



Run 


Boxsize 




muM 




e 










h-'^MQ 




D4 


33.75 


2 X 216^ 


2.75 X 10* 


4.24 X 10^ 


6.25 


D6 


33.75 


2 X 486^ 


2.75 X lO'^ 


4.24 X 10'' 


2.73 


E6 


50 


2 X 486^ 


7.85 X lO'^ 


1.21 X lO'' 


4.12 



A'p: Total number of particles 

mDM- Mass of dark matter particles 

mgas : Initial mass of gas particles 

e: Comoving gravitational softening length 



used 

2 



The accretion rate of each black hole 
to compute the bolome tric luminosity, L — t^Mbhc 
ijShakura fc Sunvaevll973l l. Here rj is the radiative efficiency, 
and it is fixed at 0.1 throughout the simulation and this 
analysis. Some coupling between the liberated luminosity 
and the surrounding gas is expected: in the simulation 5 per 
cent of the luminosity is (isotropically) deposited as ther- 
mal energy in the l ocal black hole k ernel, acting as a form 
of feedback energy (jPi Matteo et aLlUoOSil . 

Note that to derive luminosities in specific wavebands 
(consistent with the observational constraints), we need to 
apply a bolometric correction to our Cjuasar luminosi ties. We 
apply the bolome tric correction from Hop kins et al] l|2007bl ) 
(consistent with l|Marconi et al.|[2004 ')): 



Ci 



+ C2 



(1) 



where cl =(6.25, 17.87, 10.83), c2 = (9.00, 10.03, 6.08), 
kl = (-0.37, 0.28, 0.28), k2 = (-0.012, -0.020, -0.020) for B- 
Band, 0.5-2 keV Soft X-ray band, and 2-10 keV Hard X-ray 
band, respectively. 




Log(M/Me) 



Figure 1. Relation between mass and bolometric luminosity for 
black holes in the D6 simulation for redshifts 1, 3, and 5. The 
lines show LEdd (solid pink) and O.OlLgdd (dashed green). 



2.2 Simulation parameters 

Three simulation runs are analysed in this paper to allow 
testing for resolution effects. The main parameters are listed 
in Table [T] The three runs were of moderate volume, with 
boxsizes of side length 33.75/i~^Mpc (D6 and D4 simula- 
tions), and 50/i"^Mpc (E6). For the D6 and E6 runs Np = 

2 X 486^ particles were used, and the D4 used 2 x 216^. The 
moderate boxsizes prevent the simulations from being run 
below 2: ~ 1 (z 0.5 for D4 run) to keep the fundamental 
mode linear, but provide a large enough scale to produce 
sufficiently luminous sources, albeit rare. The limitation on 
the boxsizes is necessary to allow for appropriate resolution 
to carry out the subgrid physics in a converged regime (for 
further details on the simulation methods, pa rameters and 
convergence studies see lPi Matteo et al.l (j2008l ) and also the 
discussion at the end of this paper). 

3 RESULTS 

3.1 Mass-Luminosity Relation 

In order to illustrate the range of properties of the black hole 
population in the simulations, in Fig. [1] we show the relation 
between black hole mass and luminosity for the whole sam- 
ple of objects in the D6 simulation (at 2; = 1, 3, and 5). Note 




12 3 4 5 6 7 



Figure 2. Example of lightcurves between z=l and z=7 for three 
massive black holes in the D6 simulation. The top panel shows 
the growth of the black hole masses and bottom three panels (in 
corresponding colors) show the associated lightcurves for their 
bolometric luminosity. 
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the D4 and E6 mass-luminosity relations are not plotted, but 
both produced similar results. There is some correlation, al- 
beit weak, between luminosity and mass of black holes, how- 
ever in most regimes a significant scatter is seen, implying 
a fair range of luminosities for a fixed black hole mass. This 
is the direct result of our simulations and in particular the 
complex lightcurves associated with the accretion history 
(and the evolution of the gas supply) which is followed in 
detail for all the black holes in our simulations. As an exam- 
ple, in Fig. [2]we show the black hole mass assembly history 
for three specific black holes in the D6 simulations and their 
associated lightcurves (in terms of bolometric luminosity). 
The high level of variability in the lightcurves is induced by 
the detailed hydrodynamics, interplay between gas inflows, 
associated accretion and feedback p rocesses self-consi s tently 
modelled in the simulations (see also lDi Matteo et al.l (|2008l ) 
for more examples of accretion and merger histories of black 
holes in the simulations). This implies that in turn we ex- 
pect the same black holes, at different stages of activity, to 
contribute to different regions of the luminosity function. 



_ ^ p-. — , — . — p-. — . — , — I — .— . — . — p-. — . — ' I ' ' — ' — I — ' — >— ' — r 




5 8 10 12 , 14 , 8 10 12 14 

Log(L/Lo) 

Figure 3. The luminosity function as computed for the D6 sim- 
ulation at redshifts 1 and 3 for black holes restricted to the 
following ranges: Blue: < IO'^Mq; Green: 10® - 10^ Mq; Red: 
10^ - 10^M(7); Brown : > lO^M©; Black: Full mass range. The 
iHopk ins et al. best-fitting bolometric QLF has also been 

plotted for comparison (Solid black line). 



3.2 Luminosity Functions 

To illustrate the effect of different black hole populations 
to the QLF in detail, in Figure [3] we plot the relative con- 
tribution from the different black hole mass ranges to the 
luminosity function at z — 1 and 3. At z = 1, the black 
holes with mass below lO^M© provide the dominant con- 
tribution to the luminosity function for luminosities below 
lO''-^ 1/0 (10"*^ ergs" ^), while higher masses dominate at larger 
luminosities. At higher redshift, the low mass black holes 
are the dominant contribution up to a higher luminosity 
(1O^°-^L0 at 2 = 3), and have a more significant contribu- 
tion to higher luminosities than they do at low redshift. Note 
that black holes with masses < IO'^Mq give rise to a signifi- 
cant steepening of the luminosity function below ^ 10^" Lq 
(although this is sometimes below the current observational 
limits, it is an important effect in our results). In our nu- 
merical simulations we do not expect to resolve the accretion 
history on to the lowest mass black holes (as also shown by 
increased scatter in Fig. [T]) , where the gas dynamics are well 
resolved only well beyond the black hole accretion radius. 
For this reason, as well as the fact that the low-mass black 
holes correspond primarily to recently inserted seed parti- 
cles which have yet to undergo critical accretion phases (i.e. 
dependent on our initial choice of this parameter), we will 
use only black holes with Mbh > 10® Mq for the rest of our 
analysis. 

Fig. |4] shows the predictions from our simulations for 
the AGN luminosity functions for z = 0.5, 1, 2, 3, 4, 5, (note 
that only the D4 simulation is used at z — 0.5). The first 
column shows the bolometric luminosity function derived di- 
rectly from the simulations. The second and third columns 
show the luminosity function after applying the bolomet- 
ric correction (Eq.l) to obtain the B-band (second column) 
and 2-10 keV hard X-ray band (third column) QLFs. Along 
with our predictions, we plot the observational data and the 
best-fitting QLFs from several studies for the hard X-ray 



and optical bands (see Fig. |4] caption for complete list of 
references 

The bolometric luminositi es are compare d to the best- 
fitting function computed by I Hopkins et al. I (l2007b!) who 
compiled all available data from observational studies across 
several bands, including the optical, mid-infrared, hard and 
soft X-ray ba nds, and fitted them t o a double power law 
function (see iHopkins et al.l (l2007lj ) for the function and 
the table of redshift-dependent parameters). The best-fitting 
function is plotted consistent with the range of observed 
bolometric luminosities. This is why the minimum luminos- 
ity shown in these fit functions is redshift-dependent, and 
ranges from 9.97Lq at z=l to 11.53L0 at z—6. 

Comparing observed and predicted LFs, it is apparent 
that the simulations can only reproduce the 'faint-end' of 
the LF: this is expected as the number density of AGNs 
in the 'bright-end' is simply too low to be accessible in our 
simulated volumes. Thus our predictions are limited to a rel- 
atively small range of luminosities which can be compared 
directly to observational data, and the largest overlap is in 
the X-ray band, rather than in the B-band due to the signif- 
icantly fainter AGN populations in the former. Related to 
this is the lack of predictions for the knee of the QLF, which 
occurs at a higher luminosity than the simulations produce. 

Within the range of comparison, our predictions agree 
well with the data. Our simulations are fully consistent with 
the constraints from the B-band (albeit with very limited 
region of overlap). In the hard X-ray band, at L ~ 1O^''L0 
(L ~ 10*^'^ erg s~^), close to the maximum luminosities 
probed with our simulations, there is also very good agree- 



1 Note that both lEbrero et al] | |2009|) and lYencho et all (|2009|) 
considered 2-8 keV rather than 2-10 keV, so their functions were 
adjusted using a photon index of F = 1.8 to maintain a con- 
sistent defini t ion of the hard X-ray band. In addition, neither 
lEbrero et al] | |2009|) nor lYencho et al. I ll200St) cons i der a bsor] )- 
tion, whereas Ueda et al. l|2003l ). iLa Franca et al] (120051 ). and 
[Silverman et al. I ll2008h all use absorption corrected data. 
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Log(L) [erg/s] 

42 43 44 45 46 47 42 43 44 45 46 47 42 43 44 45 46 47 




8 10 12 14 10 12 14 10 12 14 

Bolometric B — Band Hard X — Ray 

Log(L/LQ) 

Figure 4. The black hole luminosity function for all three simulations (blue - D6 simulation; red - D4 simulation; green - E6 simulation) 
using sources with Mbh > IO^Mq. The first c olumn is the bolometri c QLF computed directly from the simulation. The solid black line 
is the double-power law QLF function given in lHopkins et all (l2007bl) . The second and third columns show the luminosity function after 
applying a bolometric correction (Eq. [T]l t o produce the B-Band and H ard X-ray band . Open circles for opti cal bands are dat apoints 
from the f ollowing studies: Bright purple - iRichards et all (l2005h: Blue - 'Richard s et al.l l|2006h : D ark Green - IWolf et al.l l|2003l '): Red - 
iHunt et al . (2004); Yellow - Cristia ni et al.l ll2004f): Orange -iKennefick e t al. (199j); Dark Purple - ISchmidt et al.l \l99^ : Bright Green 
- iFan et al.l (|2001a b, 2003, 2004); Black -'Sianaet al.' ^2000)7 Closed cir cles for hard X-ray bands are datapoi nts from the followin g 
studies: Pink - jUeda et al.. (.2003,) : Blue - Silverman et al. (200^!); Green - iBarger e t al. (2005, 2003.a Ji); O range - iNandra et all ( l2005h . 
Dotted lines in the har d X-ray c olumn are best-fit ting LDDE function s from th e following st udies: Pink - | Ueda et al.l ( 20031: Orange - 
iLa Franca et al.l ||2005|) ; Purple - ISilverman et al.l [2008,') : Black - .Ebrero et al.l ( liooa); Red -lYencho et al.l ||2009| ) . Dotted line s in th e 
B-Band column are be s t-fitti ng LDDE functions from the following studies: Pink - iBovle etal .1 l|2000l ): Orange - ICroom et al ] (|2004h : 
Green - IRichards et al.l (|2005[) . The dashed line in the hard X-ray at z=2 is the D6 QLF if only A^bh > IO'-^Mq are included. 
© 200? RAS, MNRAS OOO.llHIOl 
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ment. For z ^ 1 the overall shape and slope of the faint-end 
is also reproduced remarkably well. At 2 2 ^ 4 however, 
the slope predicted from the simulation is typically steeper 
than in the observed LFs. At z > 4, if compared to the 
fits of the observations, the slopes are again consistent. The 
same result is found in the comparison with the bolomet- 
ric luminosity functions (where indeed the hard X-ray data 
significantly dominates Hopkins' fits in the low luminosity 
end) and where the discrepancy in the slope is the greatest 
at 2: ~ 2. 

It is promising that the simulations agree well with the 
data at 2 ^ 1 where the observations in the 'faint-end' are 
most complete and the bolometric corrections (which are 
derived from the local samples and have no redshift depen- 
dence) are most appropriately applied. The larger number of 
AGNs at the faintest luminosities (hence the steeper slopes) 
at higher redshifts may have two possible implications. One 
is that there is still a population of the faint, possibly heav- 
ily obscured AGN above z = 2 that have not yet been de- 
tected. Alternatively, our simulations are actually overpro- 
ducing the faint AGN population due to lack of appropri- 
ate resolution or appropriate feedback physics (either due 
to stellar processes or AGN) in the faint, low mass black 
hole population. The latter possibility is also hinted at from 
our results in Fig. [3] and the fact that the D6 predictions 
(highest resolution simulation) typically predict the fiattest 
slopes (albeit barely as the convergence between the predic- 
tions for the three simulations is good). We will investigate 
this further in the course of the paper. 

3.2.1 Model dependent effects on the predictions for QLF 

Even though our simulations predict directly an accretion 
rate for each of the BHs (from the hydrodynamics) , at any 
given time there are two major model dependent assump- 
tions that we have made in order to the translate the ac- 
cretion rate into a luminosity in a particular waveband: one 
is the assumption of a fixed radiative efficiency of ten per 
cent, and the second is the use of an empirically derived 
bolometric correction. Here we wish to investigate the main 
effects on the QLF predictions of varying these two assump- 
tions. In Figure [S] we demonstrate these effect by showing 
the B-Band QLF at z=l (left) and the hard X-ray band 
QLF at z=2 (right). With these two panels we are able to 
fully illustrate the relevance of the effects. 

When calculating black hole luminosities (as described 
in Section 2.1), we are assuming that all sources have 
equal radiative efficiency and set it to 0.1. However, it has 
been suggested (even though the precise changes in accre- 
tion physics at low accretion rate remain somewhat uncer- 
tain) that sources accreting at a sufficiently low Eddington 
fraction (typically at ~ 0.01 Eddington) are expected to 
transition to a radiatively inefficient state with associated 
changes in the spectral energy distribution (SED) that are 
dramatically diff erent. These radiatively inefficient accre- 
tion models (e.g. 



lYuan fc NaravanI 



Naravanllioosi : iQuataert fc Naravan|[l999l : 
2004 and references therein) (and also ob- 



servations of both AGN and X-ray transients) indicate that 
such transitions occur around accretion rates of 1% the Ed- 
dington value. The simplest way to investigate the overall 
effect this may have on our QLF predictions is to eliminate 
all sources accreting below O.OlLEdd (blue line). As shown 



in Figure [T] most sources are above this cut-off luminos- 
ity at 2 > 1, so we expect a minimal effect on the QLF 
above this redshift. Eliminating low luminosity sources for 
2^1 leads to a flattening of the QLF slope (Fig. [Sjl, as 
most sources are actually below this threshold at this point 
(as indeed generally expected that such modes of accretion 
will be well below the quasar peak). In short, low radiative 
efficiency accretion is expected to lead to some flattening 
of the QLF at 2 ^ 1. This effect therefore would not help 
flattening the QLF function at higher redshift and therefore 
better reconcile our predictions with observations. 

Our predictions for the various wavebands are of course 
dependent on the form of the bolometric correction used (the 
one adopted here is shown in Eq. 1). Even though there is 
a luminosity dependence in our bolometric correction, it is 
less well-constrained at low luminosities (also for the rea- 
sons discussed above), where the majority of our sources lie. 
To explore the effects that the luminosity dependence has 
on our results in Figure [S] (green line) we show the QLFs 
for a luminosity-independent bolometric correction, using 
the value of the correction factor evaluated at L = 10^^ Lq, 
where the correction factor is best constrained. Doing this 
has a small effect on the B-Band QLF, where in any case (see 
Eq.l) the correction has a small dependence on luminosity. 
In the hard X-ray band, however, a luminosity independent 
correction produces significantly lower magnitude for the 
QLF, which more closely matches the observational data. 
This illustrates that the exact form of the bolometric cor- 
rection, particularly the form of its luminosity dependence, 
may have a strong effect on our final results. 



3.3 Comoving Number Density Evolution 

The quasar comoving number density evolution is plotted 
in Fig. [S] This is derived by integrating the luminosity 
functions plotted in Fig. |4] (we average over all three sim- 
ulations). Again we plot the predictions for the bolomet- 
ric, B-band, hard X-ray and in this case soft X-ray band 
also. For the latter two bands we a l so sh ow the obser- 
vational constraints from lUeda et all (|2003D (hard X-ray) 
and[lla singcr ct al. (200^) (soft X-ray). Note that, following 
I Hopkins et al.l (|2006bf ). the normalization of the soft x-ray 
data has been multiplied by 10 to adjust for obscuration in 
this band (this adjustment is somewhat model dependent. 



^ Additionally, because our simulations use a single feedback 
model for all black holes, they do not model separate 'quasar' 
and 'radio' modes. In addition to having an effect on the radia- 
tive efficiency used to determine the BH luminosity, the inclusion 
of a radio mode will have a quenching effect during the simu- 
l ation, due to the r adio mode suppressing inflow of cooling gas 
( Croton et al.|[2006h . As the majority of our sources at 2 > 1 
are accreting above 0.01 times the Eddington accretion rate, it 
is only at low redshifts (at or below 2 1, which is the limit of 
our simulations) that we would expect the radio m ode to have a 
signifi cant effect on black hole growth. Additionally. [Siiacki et al] 
ll2007h found that, although the effect of the radio mode does be- 
come large at 2 < 1 , the bulk of black hole growth is always during 
the quasar mode (with the quasar mode accretion contributing 
95% of the integrated black hole mass density) , and that model- 
ing a separate radio mode has negligible effect on Mbh — and 
A/bh — relations. 
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Figure 5. The QLF for D6 simulation for B-Band at z=l (left) and Hard X-ray at z=2 (right) using different model parameters: Black 
- same parameters as Fig. |4] Blue - sources with L > O.OlLgtjd; Green - luminosity-independent bolometric correction. See Fig. |4] for 
references of observational datapoints and best-fitting functions. 




Figure 6. The comoving number density evolution of quasars for Mbh > 10^ Mq averaged across the three simulations (except at 
2 = 0.5, where only the D4 simulation was available). A linear extrapolation has been applied to the simulation to a llow our lum i nosity 
function to be integrated over the given luminosity rang es, and be compared directly to observational data from lUeda et al.l ||2003| ) 
(hard X-ray - purple trian gles) andlHasinger et al.l l|2005l ) (soft X-ray - green tria ngles). We have al so plotted the predictions from the 
best-fitting functions from [Richards et al.l |200^ (optical - green dashed line) and lUeda et al.l l|2003l ) (hard X-ray - purple dashed line) . 
Note the optical observational data were limited to redshifts below 2 = 3, so we ha ve only plotted those curves for 2 < 3. Additionally, 
the soft X-ray data have been adjusted for obscuration as in iHopk ins et al ] l|2006bh . 



but provides a first approximation for comparison) . We have 
al so plotted the p r edicti ons from the be st-fitting fu n ction s 
of iRichards etlo] (j2005h (B-Band) and lUeda et al.l (|2003l ) 
(Hard X-ray). The B-Band function is terminated at 2; = 3 
since the fits were based only on sources below this redshift. 
In some cases a linear extrapolation (to higher luminosities) 
was applied to the simulation to allow the range of inte- 
gration to match the observational data (given in specific 
luminosity bins). 

In virtually every band, the quasar number density from 
the simulations peaks at z ~ 2.5 and as expected, their 
number density is dominated by the lower luminosity pop- 



ulations. When comparing to the X-ray data (both hard 
and soft bands), we again find there is good agreement with 
the observed evolution in the intermediate/high luminosity 
ranges. Consistent with the results from the luminosity func- 
tions, the evolution in the lowest luminosity range implies 
larger number densities in t he soft X-ra y band above 2. 
The hard X-ray data from lUeda et all l|2003l ) for the low- 
est luminosity range is limited to 2; < 0.5, preventing direct 
comparison, however the best-fitting function's extrapola- 
tion shows that we may have a similar overestimate for the 
hard X-ray band. 
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Figure 7. The total luminosity density evolution averaged across 
the three simulations (except at z = 0.5 where only the D4 sim- 
ulation was used). Top panel - bolometric prediction computed 
directly from simulation; remaining plots - optical and hard X- 
ray predictions obtained by applying bolometric corrections to 
the bolometric luminosities. Grey crosses show the luminosity 
density evolution below a certain cutoff luminosity (chosen at 
approximately the lowest break luminosity for the given band to 
restrict our predictions to the faint-end). The remaining (colored) 
symbols show the emissivity produced by black holes within the 
following luminosity ranges: Blue diamonds - 10''^ ergs~^ < L < 
lO*'^ erg s~^; Green triangles - 10^* erg s~^ < L < 10*^ erg s~^; 
Red squares - 10*'^ erg s~^ < L < 10** erg s~^; Purple X - 
10*2 erg s-i < L < 10*^ erg s-^ The shaded regions are the ar- 
eas bo u nded by the best - fitting LD DE functions f r omlUeda et al] 
ll2003l '). lLa"Franca et all ||2005|) andlEbrero et al.l ||2009|) for hard 
X-ray and for optical, [Richards et al.l | |2005^ and iBovle et all 
1I2OOOI '). Note the observational estimates for the optical band are 
only plotted for 2: < 3 since the observational data was limited to 
low redshifts. 



3.4 Luminosity Density Evolution 

In Fig. [3 we show the total luminosity density evolution 
from the simulated quasars with the appropriate observa- 
tional constraints. A linear extrapolation of the QLF was 
made for the highest luminosity bins such that consistent lu- 
minosity ranges could be used across all simulations and red- 
shifts. To cover the full range of observational constraints, 
we have used the best-fitting luminosity-dependent density 
evolution (LDDE) functions from th e following studies to 
bou nd the shaded regio ns: in optical - [Richards et al.l (|2005| ) 
andlBovle et al.l llgOOO!); in hard X-ray -lUeda et al] (12003 



iLa Franca et al.l ^05) and lEbrero et al.l (|2009l ). For these 
reasons. Fig. [7] is less of a direct comparison between simu- 



Figure 8. The total bolometric luminosity density evolution av- 
eraged across the three simulations (except at z = 0.5 where only 
the D4 simulation was run), broken into mass bins. The upper 
plot includes all black holes, while the lower plot has neglected 
the outlying black hole in the D6 simulation at redshift 3. 



lations (where we extrapolate somewhat) and observations 
(where we used model dependent fits to data as constraints) . 

While the predicted peak in the total luminosity den- 
sity is at 2 ~ 2.5 across the various bands, there are some 
fairly marked differences in the objects that dominate the 
contribution to the total luminosity evolution. The bolo- 
metric luminosity density (top panel) is dominated by the 
highest luminosity bins, with the brightest objects peaking 
at z ~ 4 while the lower luminosity bins peak at z ~ 2. In 
the optical band (second panel), the luminosity density is 
still dominated by the brightest objects, comparable to the 
observational constraints. 

The most significant difference is in the X-ray band, 
where the low luminosity population produces most of the 
contribution to the luminosity density. Additionally, the 
peak in observed luminosity density is close to z ~ 1 rather 
than 2 ~ 2.5 as implied by the simulations. The overall re- 
versal of trends in the relative contributions in various bands 
is of course caused by the form of the bolometric correction 
used in Eq.[T] 

Finally, Fig. [8] shows the contribution to the bolomet- 
ric luminosity density evolution as a function of black hole 
mass. The lower plot in Fig. [8] is identical to the upper plot 
except the single most luminous black hole from the D6 sim- 
ulation at redshift 3 has been neglected (to explicitly show 
effects due to small statistics). We find that it is typically 
the midrange black holes masses (10*^ Af© < A4"bh < IQ^Mq) 
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which provide the largest contribution to the luminosity den- 
sity. As one might expect, at higher redshifts, the low mass 
black holes provide a more significant contribution, which 
was expected due to the lack of high mass black holes (see 
also Fig. [T|. 



4 DISCUSSION 

One of our primary results from the luminosity function and 
the global number density and luminosity density evolution 
is that our simulations are in good agreement with the ob- 
servational constraints but imply a larger number of low 
luminosity X-ray sources (at 2 > 2) than observed. In or- 
der to assess the viability of our results we need to check 
whether this population may violate the current constraints 
on the total X-ray backgr ound (XRB). T he XRB intensity 
is calculated according to ()Peacocklll999l ) 

where the emissivity, is the hard X-ray emissivity 
shown in Fig. [T] z is the redshift, [21,22] is the range of 
redshifts being considered, z/q is the frequency at redshift 
zero, Hq is the Hubble Parameter at redshift zero, and f2o 
is the total density parameter at redshift zero. A photon 
index of F = 1.8 was assumed when computing e^([l -I- z\vq) 
to account for the form of the power spectrum of the black 
holes. Since we only have simulation information at discrete 
redshifts, we interpolate linearly between datapoints to 
compute the integral. 

We find that the total contribution to the 2-10 
keV X-ray background from our simulated black holes is 
/2-iokev,D6 = 1-28 X 10""ergs"^cm"^deg"^ for the D6 
simulation (recall the simulations are restricted to 2 ^ 
1). This is well within the observed 2-10 keV hard X- 
ray background intensit y, J2-iokeV.obs = 2 .02 ±0.11 x 
10""ergs"^cm"^deg"^ (|Moretti et al.l 120031 '). If we apply 
a linear extrapolation to the simulation to include an ap- 
proximate contribution from 2 < 1, we predict a total XRB 
intensity of /2-iokcV = 1-8 x lO^^^erg s^^cm^^deg"^, still 
below the observed value. 

When we do a similar calculation using the E6 and 
D4 simulations (which have lower resolution) we produce 
XRB intensities of /2-iokoV = 1-94 x 10"" and 2.87 x 
10~^^erg s~^cm~^deg~'^ for the E6 and D4 simultions at 
2 > 1, respectively. With the D4 simulation we are starting 
to violate the observed XRB even without including 2 < 1. 
This result is quite fundamental for our analysis as it indi- 
cates that the simulation resolution plays an important role 
in estimating the effect from the faintest sources (those that 
were found to be in excess of the observed LF) and that 
their contribution decreases with higher resolution (note in 
Fig. m the steepest slopes of LF are always predicted from 
the D4 run). 

Our simulations also predict that the QLFs extend to 
fainter fiuxes than currently observed. We therefore com- 
pare our predictions for the unresolved fraction of the 2-8 
keV XRB from the D6 simulation to test whether current 
observational constraints are still consistent with our predic- 
tions (i.e. if we could still be missing a faint population of 
e.g. heavily obscured AGN). We found that our prediction 



for the XRB contribution from luminosities limited to the 
range of overlap between observation and simulation pro- 
vides an excess intensity of 1.85 x 10~^^erg s~^cm~^deg~^ 
in the 2-8 keV band (assuming photon index of 1.8). This is 
below the 2-8 keV unresolved backgro und of 72-8keV.unres = 
3.4 ± 1.7 X 10"^^ergs"^cm"^deg"^ l|Hickox fc MarkevitchI 
l2006h . If we take into account the total contribution from 
the whole population in the simulations (well below the 
faintest sources currently observed but still assuming the 
same X-ray spectrum) we are in excess of the unresolved 
background by almost a factor of 2. To further illustrate 
and elucidate this issue, in Fig. [5] we plot the differential 
contribution to the 2-10 keV background from the simu- 
lations and the observations in several redshift bins (the 
filled curves are the regions bounded by th e pre dictions from 
lUeda et al.l (120031 ') . iLa Franca et al.1 (|2005l ) and lEbrero et all 
(20oJ)7'This shows again that the excess in our predictions 
is caused by the contribution from low-luminosity sources 
(L < lO'^'^erg s~^) and originates mostly at 2 ^ 2. This sup- 
ports the idea that the faintest sources at high redshift are 
problematic in our predictions. Figs. 3 and 7 do indeed show 
that above z=2 this population is dominated by the lowest 
mass black holes (as opposed to lower redshifts where a more 
significant fraction of high mass black holes have 'turned- 
off') and those that are likely to suffer more strongly from 
lack of resolution. For illustration, in Fig.|4l at z = 2, in the 
hard X-ray band we show how the predictions look when 
only A/bh > IO'^ '^Mq are plotted. The excess in our predic- 
tion is eliminated and the lowest luminosity end of the LF 
is now in good agreement with all the observations. 

Note also that we have used a redshift-independent cor- 
rection to convert the simulations' bolometric luminosity to 
the hard X-ray band to compare with observations. Direct 
determinations of bolometric the correction as a function of 
redshift are not yet available and this may further bias our 
results at 2 ^ 2. 



5 CONCLUSIONS 

Here we study the luminosity function and its evolution for 
populations of quasars extracted from full cosmological hy- 
drodynamical simulations which include direct modelling for 
the growth of black holes. Noting that our simulations (due 
to limitations on the volumes probed) can only be used to 
study the faint-end of the QLF, we summarize our main 
results as follows: 

• Consistent with the complex light-curves and various 
phases of activity that black holes undergo through their 
cosmic history, we have shown that there is a significant 
spread in luminosities for a given black hole mass, and in 
turn that different black hole masses contribute to the same 
regions of the QLF. 

• At low redshift (2 < 2), the low luminosity ranges (be- 
low 10^ L0) of the QLF are dominated by black holes below 
lO^M©, while the luminosities above 10® contain compa- 
rable contributions from both low and high masses. At high 
redshift the majority of black holes are below lO^M©, and 
thus the entire QLF is dominated by low black hole mass 
sources. 

• We have shown that our predictions for the faint-end of 
the QLF agree remarkably well with observations at 2 ^ 1, 
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ray background, it appears that our simulations are overpro- 
ducing low luminosity sources, particularly at intermediate 
redshifts. We have shown however that the higher resolution 
simulations produce fewer low-luminosity black holes, which 
makes it likely that this overproduction is dominated by res- 
olution effects. Additionally, our results are most accurate 
at low redshift, when the high mass (and thus least likely to 
be affected by resolution limits) sources are most important, 
further suggesting that our overproduction of low luminosity 
sources is dominated by resolution effects. 

Overall our results support the interpretation of the faint- 
end luminosity function put forward by Hopkins et al. In 
upcoming work we will compare detailed characteristics of 
black hole lightcurves in our simulations and compare their 
instantaneous luminosities to their peak luminosities, so as 
to determine more precisely if the faint-end slope is domi- 
nated by quasars radiating below their peak, or by quasars 
with faint peak luminosities, as previous models assumed. It 
may be possible (although currently infeasible due to tech- 
nological constraints) to run larger volume simulations at 
similar or higher resolution to increase the statistics at the 
bright-end of the luminosity function and further investi- 
gate the rapid dropoff in comoving number density at 2 < 1 
found in the observational data. 



Figure 9. The differential contribution to the 2-10 kcV hard 
X-ray background intensity as a function of luminosity, plotted 
for several redshift bins. The shaded areas s how the region of 
contributions bounded by t he predictions by|yeda et al., (20030, 
iLa Franca et al.l bOOSi ') and lEbrero et al.l 1 I2OO9I) . The dotted line 
is an extrapolation of the simulation obtained by extending the 
luminosity density evolution from z = I to z = 0. 



but produce steeper slopes than implied by current con- 
straints for the hard X-ray band at redshifts z = 2 and 
3. 

• Taking into account a possible transition to low radia- 
tive efficiency accretion modes for low accretion rate sources 
tends to ffatten the QLF at low redshifts but this does not 
affect significantly any of our results. 

• The exact form of the bolometric correction has a sig- 
nificant effect on our predictions. In particular, when com- 
paring to a fixed correction, the empirically determined (Eq. 
1) luminosity dependence leads to a larger QLF magnitude 
in the X-ray band. Note also that in addition, no constraints 
are currently available on the redshift dependence of the cor- 
rection. 

• The evolution of the comoving number density is 
in agreement with current constraints for the luminosity 
ranges above lO^^ergs"^. Agreement for luminosities be- 
low lO^^ergs"^ is significantly worse, but the more limited 
observational data at these ranges combined with the domi- 
nance of black holes with Mbh < IO^Mq , which are less- well 
resolved in our simulation makes these results less meaning- 
ful. 

• The luminosity density evolution predicts a peak lumi- 
nosity density at z — 2.5, with comparable contributions 
from different luminosity bins. 

• Based on the slope of the faint-end QLF, the luminosity 
density evolution and a moderate excess in the unresolved X- 
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